Hands-on Exercise 4: 2nd Order Spatial Point Patterns Analysis Methods

Author

Lorielle Malveda

1. OVERVIEW

This is the 1st of 2 Hands-on Exercise devoted to Spatial Point Pattern Analysis.

Spatial Point Pattern Analysis is the evalutation of any pattern or distribution of a set of points on a surface. These points can represent locations of various events like crimes, traffic accidents, and disease outbreaks, or they could indicate business services such as coffee shops and fast food outlets, as well as facilities like childcare and eldercare centers.

In this Hands-on Exercise, we will use relevant functions from the spatstat (for more info, click here) package to explore the spatial distribution of childcare centers in Singapore.

Specifically, this exercise aims to answer the following questions:

  • Are childcare centers in Singapore randomly distributed across the country?

  • If they are not, then where are the areas with a higher concentration of childcare centers?

2. GETTING STARTED

2.1 Acquire Data

The following datasets will be used:

  • CHILDCARE in geojson format (from data.gov.sg)

    • a point feature data providing both location and attribute information of childcare centers.
  • MP14_SUBZONE_WEB_PL in ESRI shapefile format (from data.gov.sg)

    • a polygon feature data providing information of the URA 2014 Master Plan Planning Subzone boundary data.
  • CostalOutline in ESRI shapefile format (from SLA)

    • a polygon feature data showing the national boundary of Singapore.

2.2 Install and Launch R Packages

pacman::p_load(sf, raster, spatstat, tmap, tidyverse)

This code chunk uses p_load() of the pacman package (stands for Package Manager) to check if the following packages are installed:

  • sf: designed to import, manage, and process vector-based geospatial data in R.

  • spatstat: useful for point pattern analysis. We will use this to perform 1st- and 2nd-order Spatial Point Patterns Analysis and to derive the Kernel Density Estimation (KDE) layer.

  • raster: reads, writes, manipulates, analyzes, and model gridded spatial data (i.e. raster). We will use this to convert image outputs generated by spatstat into raster format.

  • maptools: provides tools for manipulating geographic data. We will mainly use this to convert Spatial objects into the ppp format required by spatstat.

  • tmap: you must now be familiar with this as we have used this in the previous Hands-on Exercises. This provides functions for cartographic quality static point patterns or interactive maps by using the Leaflet API.

If available, the packages will then be launched into R.

3. SPATIAL DATA WRANGLING

3.1 Importing the Spatial Data

To import the datasets, st_read() of the sf package will be used.

childcare_sf <- st_read("C:/loriellemalveda/ISSS626-GAA/Hands-on_Ex/Hands-on_Ex04/data/data/child-care-services-geojson.geojson") %>%
  st_transform(crs = 3414)
Reading layer `child-care-services-geojson' from data source 
  `C:\loriellemalveda\ISSS626-GAA\Hands-on_Ex\Hands-on_Ex04\data\data\child-care-services-geojson.geojson' 
  using driver `GeoJSON'
Simple feature collection with 1545 features and 2 fields
Geometry type: POINT
Dimension:     XYZ
Bounding box:  xmin: 103.6824 ymin: 1.248403 xmax: 103.9897 ymax: 1.462134
z_range:       zmin: 0 zmax: 0
Geodetic CRS:  WGS 84
sg_sf <- st_read(dsn = "C:/loriellemalveda/ISSS626-GAA/Hands-on_Ex/Hands-on_Ex04/data/data", layer="CostalOutline")
Reading layer `CostalOutline' from data source 
  `C:\loriellemalveda\ISSS626-GAA\Hands-on_Ex\Hands-on_Ex04\data\data' 
  using driver `ESRI Shapefile'
Simple feature collection with 60 features and 4 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 2663.926 ymin: 16357.98 xmax: 56047.79 ymax: 50244.03
Projected CRS: SVY21
mpsz_sf <- st_read(dsn = "C:/loriellemalveda/ISSS626-GAA/Hands-on_Ex/Hands-on_Ex04/data/data", 
                layer = "MP14_SUBZONE_WEB_PL")
Reading layer `MP14_SUBZONE_WEB_PL' from data source 
  `C:\loriellemalveda\ISSS626-GAA\Hands-on_Ex\Hands-on_Ex04\data\data' 
  using driver `ESRI Shapefile'
Simple feature collection with 323 features and 15 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 2667.538 ymin: 15748.72 xmax: 56396.44 ymax: 50256.33
Projected CRS: SVY21

Before proceeding to analysis, let us check first if they are all in the same projection system. Do not forget to do this! :)

REVIEW: Using the sf function we learned before, we will retrieve the referencing system information of the geospatial data.

st_crs(childcare_sf)
Coordinate Reference System:
  User input: EPSG:3414 
  wkt:
PROJCRS["SVY21 / Singapore TM",
    BASEGEOGCRS["SVY21",
        DATUM["SVY21",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4757]],
    CONVERSION["Singapore Transverse Mercator",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",1.36666666666667,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",103.833333333333,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",1,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",28001.642,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",38744.572,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["northing (N)",north,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["easting (E)",east,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Cadastre, engineering survey, topographic mapping."],
        AREA["Singapore - onshore and offshore."],
        BBOX[1.13,103.59,1.47,104.07]],
    ID["EPSG",3414]]
st_crs(sg_sf)
Coordinate Reference System:
  User input: SVY21 
  wkt:
PROJCRS["SVY21",
    BASEGEOGCRS["SVY21[WGS84]",
        DATUM["World Geodetic System 1984",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]],
            ID["EPSG",6326]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["Degree",0.0174532925199433]]],
    CONVERSION["unnamed",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",1.36666666666667,
            ANGLEUNIT["Degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",103.833333333333,
            ANGLEUNIT["Degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",1,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",28001.642,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",38744.572,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]]]
st_crs(mpsz_sf)
Coordinate Reference System:
  User input: SVY21 
  wkt:
PROJCRS["SVY21",
    BASEGEOGCRS["SVY21[WGS84]",
        DATUM["World Geodetic System 1984",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]],
            ID["EPSG",6326]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["Degree",0.0174532925199433]]],
    CONVERSION["unnamed",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",1.36666666666667,
            ANGLEUNIT["Degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",103.833333333333,
            ANGLEUNIT["Degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",1,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",28001.642,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",38744.572,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]]]

Notice that except for childcare_sf, mpsz_sf and sg_sf do not have proper crs information.

REVIEW: Remember what was discussed in Lesson 2 and assign the correct crs to the 2 datasets.

mpsz_sf3414 <- st_transform(mpsz_sf, 
                              crs = 3414)
sg_sf3414 <- st_transform(sg_sf, 
                              crs = 3414)

Checking!

st_crs(mpsz_sf3414)
Coordinate Reference System:
  User input: EPSG:3414 
  wkt:
PROJCRS["SVY21 / Singapore TM",
    BASEGEOGCRS["SVY21",
        DATUM["SVY21",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4757]],
    CONVERSION["Singapore Transverse Mercator",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",1.36666666666667,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",103.833333333333,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",1,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",28001.642,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",38744.572,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["northing (N)",north,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["easting (E)",east,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Cadastre, engineering survey, topographic mapping."],
        AREA["Singapore - onshore and offshore."],
        BBOX[1.13,103.59,1.47,104.07]],
    ID["EPSG",3414]]
st_crs(sg_sf3414)
Coordinate Reference System:
  User input: EPSG:3414 
  wkt:
PROJCRS["SVY21 / Singapore TM",
    BASEGEOGCRS["SVY21",
        DATUM["SVY21",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4757]],
    CONVERSION["Singapore Transverse Mercator",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",1.36666666666667,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",103.833333333333,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",1,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",28001.642,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",38744.572,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["northing (N)",north,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["easting (E)",east,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Cadastre, engineering survey, topographic mapping."],
        AREA["Singapore - onshore and offshore."],
        BBOX[1.13,103.59,1.47,104.07]],
    ID["EPSG",3414]]

There is no need to change the referencing system as it is already at the Singapore National Projected Coordinate System.

3.2 Mapping the Geospatial Datasets

After verifying the coordinate reference system for each geospatial dataset, it is also helpful to create a map to visualize their spatial patterns.

Review: Use the mapping techniques learned prior to this exercise, and create a map.

tmap_mode("plot")
tmap mode set to plotting
tm_shape(mpsz_sf3414) +  
  tm_polygons(col = "lightgrey") +
  tm_shape(childcare_sf) + 
  tm_dots()

This code sets the tmap mode to “plot” for static map visualization. It then creates a map where the mpsz_sf3414 polygon layer is displayed with light grey polygons, and the childcare_sf point layer is overlaid as dots to show the locations of childcare centers.

Notice that all the geospatial layers fit within the same map extent. This shows that their referencing system and coordinate values are within a common spatial context.

Alternatively, we can also prepare a more interactive map by using the code chunk below.

tmap_mode('view')
tmap mode set to interactive viewing
tm_shape(childcare_sf)+
  tm_dots()

At the interactive mode, tmap uses Leaflet for R’s API. With an interactive map, you can freely navigate and zoom in and out of the map. It also adds to the experience of exploring!

Additionally, you can change the map’s background layer. Currently, three internet map layers are provided. They are ESRI.WorldGrayCanvas, OpenStreetMap, and ESRI.WorldTopoMap. The default is ESRI.WorldGrayCanvas.

 tmap_mode('plot')
tmap mode set to plotting

Reminder! Always remember to switch back to plot mode after the interactive map.

Interactive mode uses a connection, and you should avoid displaying too many interactive maps (no more than 10) in a single RMarkdown document when publishing on Netlify.

4. GEOSPATIAL DATA WRANGLING

Although simple feature data frames are becoming increasingly popular compared to sp’s Spatial* classes, many geospatial analysis packages still require input data in sp’s Spatial* format. In this section, we are going to convert a simple feature data frame to sp’s Spatial* class.

4.1 Converting from sf format into spatstat’s ppp format

The code chunk below uses the ppp() function of the spatstat package to convert spatial data to spatstat’s ppp object format.

childcare_ppp <- as.ppp(childcare_sf)
Warning in as.ppp.sf(childcare_sf): only first attribute column is used for
marks
childcare_ppp
Marked planar point pattern: 1545 points
marks are of storage type  'character'
window: rectangle = [11203.01, 45404.24] x [25667.6, 49300.88] units

Now, let us plot childcare_ppp and examine the difference.

plot(childcare_ppp)
Warning in default.charmap(ntypes, chars): Too many types to display every type
as a different character
Warning: Only 10 out of 1545 symbols are shown in the symbol map

Let us also take a look at the summary statistics to get an overview/feel of how the data is like.

summary(childcare_ppp)
Marked planar point pattern:  1545 points
Average intensity 1.91145e-06 points per square unit

Coordinates are given to 11 decimal places

marks are of type 'character'
Summary:
   Length     Class      Mode 
     1545 character character 

Window: rectangle = [11203.01, 45404.24] x [25667.6, 49300.88] units
                    (34200 x 23630 units)
Window area = 808287000 square units

4.2 Handling duplicated points

Here is how we check for duplicates in a ppp object:

any(duplicated(childcare_ppp))
[1] FALSE

Seems there are no duplicates. We will then use the multiplicity() function to count the number of co-incidence points. According to rdocumentation.org, the function is intended to count the number of duplicates of each element of an object.

multiplicity(childcare_ppp)
   [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
  [38] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
  [75] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [112] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [149] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [186] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [223] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [260] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [297] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [334] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [371] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [408] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [445] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [482] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [519] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [556] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [593] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [630] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [667] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [704] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [741] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [778] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [815] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [852] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [889] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [926] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [963] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[1000] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[1037] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[1074] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[1111] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[1148] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[1185] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[1222] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[1259] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[1296] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[1333] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[1370] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[1407] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[1444] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[1481] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[1518] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1

To check how many locations have more than one-point event:

sum(multiplicity(childcare_ppp) > 1)
[1] 0

The output shows that there are 128 duplicated point events. To view the locations of the duplicate point events, we will plot childcare data by using the code chunk below.

tmap_mode('view')
tmap mode set to interactive viewing
tm_shape(childcare_sf) +
  tm_dots(alpha=0.4, 
          size=0.05)

Do not forget to put tmap_mode back to ‘plot’! :)

tmap_mode('plot')
tmap mode set to plotting

How to spot the duplicate points from the map shown above?

Check the interactive map carefully. The darker shaded points are those that are overlapping.

Treating Duplicates:

Resolving duplicates can be approached in three ways.

  1. The simplest method is to delete the duplicates, though this would result in a loss of someway of valuable point events.
  2. Apply jittering, which introduces slight variations to the duplicate points, preventing them from occupying the same exact position.
  3. Make each point “unique” and then attaching the duplicate points as marks or attributes. This requires using analytical techniques that account for these marks.

The code chunk below implements the jittering approach.

childcare_ppp_jit <- rjitter(childcare_ppp, 
                             retry=TRUE, 
                             nsim=1, 
                             drop=TRUE)

Checking for the new dataset where the jittering approach has been used:

any(duplicated(childcare_ppp_jit))
[1] FALSE
multiplicity(childcare_ppp_jit)
   [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
  [38] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
  [75] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [112] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [149] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [186] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [223] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [260] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [297] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [334] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [371] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [408] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [445] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [482] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [519] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [556] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [593] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [630] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [667] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [704] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [741] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [778] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [815] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [852] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [889] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [926] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [963] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[1000] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[1037] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[1074] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[1111] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[1148] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[1185] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[1222] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[1259] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[1296] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[1333] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[1370] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[1407] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[1444] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[1481] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[1518] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
sum(multiplicity(childcare_ppp_jit) > 1)
[1] 0

4.3 Creating owin object

When analyzing spatial point patterns, it is recommended to limit the analysis to a specific geographical area, for example, the boundaries of Singapore. An object called owin under the spatstat package is specifically designed to represent this polygonal region.

The code snippet below shows how to convert the sgSpatialPolygon object into an owin object in spatstat.

sg_owin <- as.owin(sg_sf)

Use the plot() function to display the output object.

plot(sg_owin)

and summary() function of Base R.

summary(sg_owin)
Window: polygonal boundary
50 separate polygons (1 hole)
                 vertices         area relative.area
polygon 1 (hole)       30     -7081.18     -9.76e-06
polygon 2              55     82537.90      1.14e-04
polygon 3              90    415092.00      5.72e-04
polygon 4              49     16698.60      2.30e-05
polygon 5              38     24249.20      3.34e-05
polygon 6             976  23344700.00      3.22e-02
polygon 7             721   1927950.00      2.66e-03
polygon 8            1992   9992170.00      1.38e-02
polygon 9             330   1118960.00      1.54e-03
polygon 10            175    925904.00      1.28e-03
polygon 11            115    928394.00      1.28e-03
polygon 12             24      6352.39      8.76e-06
polygon 13            190    202489.00      2.79e-04
polygon 14             37     10170.50      1.40e-05
polygon 15             25     16622.70      2.29e-05
polygon 16             10      2145.07      2.96e-06
polygon 17             66     16184.10      2.23e-05
polygon 18           5195 636837000.00      8.78e-01
polygon 19             76    312332.00      4.31e-04
polygon 20            627  31891300.00      4.40e-02
polygon 21             20     32842.00      4.53e-05
polygon 22             42     55831.70      7.70e-05
polygon 23             67   1313540.00      1.81e-03
polygon 24            734   4690930.00      6.47e-03
polygon 25             16      3194.60      4.40e-06
polygon 26             15      4872.96      6.72e-06
polygon 27             15      4464.20      6.15e-06
polygon 28             14      5466.74      7.54e-06
polygon 29             37      5261.94      7.25e-06
polygon 30            111    662927.00      9.14e-04
polygon 31             69     56313.40      7.76e-05
polygon 32            143    145139.00      2.00e-04
polygon 33            397   2488210.00      3.43e-03
polygon 34             90    115991.00      1.60e-04
polygon 35             98     62682.90      8.64e-05
polygon 36            165    338736.00      4.67e-04
polygon 37            130     94046.50      1.30e-04
polygon 38             93    430642.00      5.94e-04
polygon 39             16      2010.46      2.77e-06
polygon 40            415   3253840.00      4.49e-03
polygon 41             30     10838.20      1.49e-05
polygon 42             53     34400.30      4.74e-05
polygon 43             26      8347.58      1.15e-05
polygon 44             74     58223.40      8.03e-05
polygon 45            327   2169210.00      2.99e-03
polygon 46            177    467446.00      6.44e-04
polygon 47             46    699702.00      9.65e-04
polygon 48              6     16841.00      2.32e-05
polygon 49             13     70087.30      9.66e-05
polygon 50              4      9459.63      1.30e-05
enclosing rectangle: [2663.93, 56047.79] x [16357.98, 50244.03] units
                     (53380 x 33890 units)
Window area = 725376000 square units
Fraction of frame area: 0.401

4.4 Combining point events object and owin object.

We are now at the last step of geospatial data wrangling! Check out the code below to extract childcare events that are located within Singapore.

childcareSG_ppp = childcare_ppp[sg_owin]

The output object combined both the point and polygon feature in one ppp object class as shown below.

summary(childcareSG_ppp)
Marked planar point pattern:  1545 points
Average intensity 2.129929e-06 points per square unit

Coordinates are given to 11 decimal places

marks are of type 'character'
Summary:
   Length     Class      Mode 
     1545 character character 

Window: polygonal boundary
50 separate polygons (1 hole)
                 vertices         area relative.area
polygon 1 (hole)       30     -7081.18     -9.76e-06
polygon 2              55     82537.90      1.14e-04
polygon 3              90    415092.00      5.72e-04
polygon 4              49     16698.60      2.30e-05
polygon 5              38     24249.20      3.34e-05
polygon 6             976  23344700.00      3.22e-02
polygon 7             721   1927950.00      2.66e-03
polygon 8            1992   9992170.00      1.38e-02
polygon 9             330   1118960.00      1.54e-03
polygon 10            175    925904.00      1.28e-03
polygon 11            115    928394.00      1.28e-03
polygon 12             24      6352.39      8.76e-06
polygon 13            190    202489.00      2.79e-04
polygon 14             37     10170.50      1.40e-05
polygon 15             25     16622.70      2.29e-05
polygon 16             10      2145.07      2.96e-06
polygon 17             66     16184.10      2.23e-05
polygon 18           5195 636837000.00      8.78e-01
polygon 19             76    312332.00      4.31e-04
polygon 20            627  31891300.00      4.40e-02
polygon 21             20     32842.00      4.53e-05
polygon 22             42     55831.70      7.70e-05
polygon 23             67   1313540.00      1.81e-03
polygon 24            734   4690930.00      6.47e-03
polygon 25             16      3194.60      4.40e-06
polygon 26             15      4872.96      6.72e-06
polygon 27             15      4464.20      6.15e-06
polygon 28             14      5466.74      7.54e-06
polygon 29             37      5261.94      7.25e-06
polygon 30            111    662927.00      9.14e-04
polygon 31             69     56313.40      7.76e-05
polygon 32            143    145139.00      2.00e-04
polygon 33            397   2488210.00      3.43e-03
polygon 34             90    115991.00      1.60e-04
polygon 35             98     62682.90      8.64e-05
polygon 36            165    338736.00      4.67e-04
polygon 37            130     94046.50      1.30e-04
polygon 38             93    430642.00      5.94e-04
polygon 39             16      2010.46      2.77e-06
polygon 40            415   3253840.00      4.49e-03
polygon 41             30     10838.20      1.49e-05
polygon 42             53     34400.30      4.74e-05
polygon 43             26      8347.58      1.15e-05
polygon 44             74     58223.40      8.03e-05
polygon 45            327   2169210.00      2.99e-03
polygon 46            177    467446.00      6.44e-04
polygon 47             46    699702.00      9.65e-04
polygon 48              6     16841.00      2.32e-05
polygon 49             13     70087.30      9.66e-05
polygon 50              4      9459.63      1.30e-05
enclosing rectangle: [2663.93, 56047.79] x [16357.98, 50244.03] units
                     (53380 x 33890 units)
Window area = 725376000 square units
Fraction of frame area: 0.401

Let us now plot the derived childcareSG_ppp.

plot(childcareSG_ppp)
Warning in default.charmap(ntypes, chars): Too many types to display every type
as a different character
Warning: Only 10 out of 1545 symbols are shown in the symbol map

EXTRACTING STUDY AREA

Let us extract the target planning areas.

pg <- mpsz_sf %>%
  filter(PLN_AREA_N == "PUNGGOL")
tm <- mpsz_sf %>%
  filter(PLN_AREA_N == "TAMPINES")
ck <- mpsz_sf %>%
  filter(PLN_AREA_N == "CHOA CHU KANG")
jw <- mpsz_sf %>%
  filter(PLN_AREA_N == "JURONG WEST")

Let us next plot the planning areas for better visualization

par(mfrow=c(2,2))
plot(pg, main = "Ponggol")
Warning: plotting the first 9 out of 15 attributes; use max.plot = 15 to plot
all

Visualizing each planning area:

plot(tm, main = "Tampines")
Warning: plotting the first 9 out of 15 attributes; use max.plot = 15 to plot
all

plot(ck, main = "Choa Chu Kang")
Warning: plotting the first 10 out of 15 attributes; use max.plot = 15 to plot
all

plot(jw, main = "Jurong West")
Warning: plotting the first 9 out of 15 attributes; use max.plot = 15 to plot
all

CONVERTING OF sf OBJECTS INTO owin OBJECTS

Now, we will convert these sf objects into owin objects that is required by spatstat.

pg_owin = as.owin(pg)
tm_owin = as.owin(tm)
ck_owin = as.owin(ck)
jw_owin = as.owin(jw)

COMBINING CHILDCARE POINTS AND THE STUDY AREA

Next, we will extract the childcare centers that are within the specific region for better analysis later on.

childcare_pg_ppp = childcare_ppp_jit[pg_owin]
childcare_tm_ppp = childcare_ppp_jit[tm_owin]
childcare_ck_ppp = childcare_ppp_jit[ck_owin]
childcare_jw_ppp = childcare_ppp_jit[jw_owin]

Let’s use the rescale() function to transform the unit of measurement from meter to kilometer.

childcare_pg_ppp.km = rescale(childcare_pg_ppp, 1000, "km")
childcare_tm_ppp.km = rescale(childcare_tm_ppp, 1000, "km")
childcare_ck_ppp.km = rescale(childcare_ck_ppp, 1000, "km")
childcare_jw_ppp.km = rescale(childcare_jw_ppp, 1000, "km")

Plot the study areas and locations of the childcare centers using the code chunk below.

par(mfrow=c(2,2))
plot(childcare_pg_ppp.km, main="Punggol")
Warning in default.charmap(ntypes, chars): Too many types to display every type
as a different character
Warning: Only 10 out of 61 symbols are shown in the symbol map
plot(childcare_tm_ppp.km, main="Tampines")
Warning in default.charmap(ntypes, chars): Too many types to display every type
as a different character
Warning: Only 10 out of 89 symbols are shown in the symbol map
plot(childcare_ck_ppp.km, main="Choa Chu Kang")
Warning in default.charmap(ntypes, chars): Too many types to display every type
as a different character
Warning: Only 10 out of 61 symbols are shown in the symbol map
plot(childcare_jw_ppp.km, main="Jurong West")
Warning in default.charmap(ntypes, chars): Too many types to display every type
as a different character
Warning: Only 10 out of 88 symbols are shown in the symbol map

Voila!

5. SECOND-ORDER SPATIAL PATTERNS ANALYSIS

5.1 Analyzing Spatial Point Process Using the G-Function

The G-function measures the distribution of the distances from an arbitrary event to its nearest event.

In this section, we are using:

  • Gest() of the spatstat package to compute the G-Function Estimation

  • envelope() of the spatstat package to perform Monte Carlo Simulation Tests.

CHOA CHU KANG PLANNING AREA

Computing the G-Function Estimation

The code chunk below computes the G-Function Estimation:

G_CK = Gest(childcare_ck_ppp, correction = "border")
plot(G_CK, xlim=c(0,500))

Performing Complete Spatial Randomness Test

We will conduct a hypothesis test to confirm the observed spatial patterns above.

Ho = The distribution of childcare services in Choa Chu Kang is randomly distributed.

H1= The distribution of childcare services in Choa Chu Kang is not randomly distributed.

The null hypothesis will be rejected if the p-value is smaller than alpha value of 0.001.

Monte Carlo Test:

G_CK.csr <- envelope(childcare_ck_ppp, Gest, nsim = 999)
Generating 999 simulations of CSR  ...
1, 2, 3, ......10.........20.........30.........40.........50.........60..
.......70.........80.........90.........100.........110.........120.........130
.........140.........150.........160.........170.........180.........190........
.200.........210.........220.........230.........240.........250.........260......
...270.........280.........290.........300.........310.........320.........330....
.....340.........350.........360.........370.........380.........390.........400..
.......410.........420.........430.........440.........450.........460.........470
.........480.........490.........500.........510.........520.........530........
.540.........550.........560.........570.........580.........590.........600......
...610.........620.........630.........640.........650.........660.........670....
.....680.........690.........700.........710.........720.........730.........740..
.......750.........760.........770.........780.........790.........800.........810
.........820.........830.........840.........850.........860.........870........
.880.........890.........900.........910.........920.........930.........940......
...950.........960.........970.........980.........990........
999.

Done.
plot(G_CK.csr)

TAMPINES PLANNING AREA

We are doing the same process but for the Tampines Planning Area this time.

Computing G-Function Estimation

G_tm = Gest(childcare_tm_ppp, correction = "best")
plot(G_tm)

Performing Complete Spatial Randomness Test

Ho = The distribution of childcare services in Tampines is randomly distributed.

H1= The distribution of childcare services in Tampines is not randomly distributed.

The null hypothesis will be rejected if the p-value is smaller than alpha value of 0.001.

G_tm.csr <- envelope(childcare_tm_ppp, Gest, correction = "all", nsim = 999)
Generating 999 simulations of CSR  ...
1, 2, 3, ......10.........20.........30.........40.........50.........60..
.......70.........80.........90.........100.........110.........120.........130
.........140.........150.........160.........170.........180.........190........
.200.........210.........220.........230.........240.........250.........260......
...270.........280.........290.........300.........310.........320.........330....
.....340.........350.........360.........370.........380.........390.........400..
.......410.........420.........430.........440.........450.........460.........470
.........480.........490.........500.........510.........520.........530........
.540.........550.........560.........570.........580.........590.........600......
...610.........620.........630.........640.........650.........660.........670....
.....680.........690.........700.........710.........720.........730.........740..
.......750.........760.........770.........780.........790.........800.........810
.........820.........830.........840.........850.........860.........870........
.880.........890.........900.........910.........920.........930.........940......
...950.........960.........970.........980.........990........
999.

Done.
plot(G_tm.csr)

5.2 Analyzing Spatial Point Process Using F-Function

The F-Function estimates the empty space function F(r) or its hazard rate h(r) from a point pattern in a window of arbitrary shape.

In this section, we are using:

  • Fest() of the spatstat package to compute the G-Function Estimation

  • envelope() of the spatstat package to perform Monte Carlo Simulation Tests.

CHOA CHU KANG PLANNING AREA

Computing F-Function Estimation

The code chunk below computes the F-Function Estimation:

F_CK = Fest(childcare_ck_ppp)
plot(F_CK)

Performing Complete Spatial Randomness Test

We will conduct a hypothesis test to confirm the observed spatial patterns above.

Ho = The distribution of childcare services in Choa Chu Kang is randomly distributed.

H1= The distribution of childcare services in Choa Chu Kang is not randomly distributed.

The null hypothesis will be rejected if the p-value is smaller than alpha value of 0.001.

Monte Carlo Test:

F_CK.csr <- envelope(childcare_ck_ppp, Fest, nsim = 999)
Generating 999 simulations of CSR  ...
1, 2, 3, ......10.........20.........30.........40.........50.........60..
.......70.........80.........90.........100.........110.........120.........130
.........140.........150.........160.........170.........180.........190........
.200.........210.........220.........230.........240.........250.........260......
...270.........280.........290.........300.........310.........320.........330....
.....340.........350.........360.........370.........380.........390.........400..
.......410.........420.........430.........440.........450.........460.........470
.........480.........490.........500.........510.........520.........530........
.540.........550.........560.........570.........580.........590.........600......
...610.........620.........630.........640.........650.........660.........670....
.....680.........690.........700.........710.........720.........730.........740..
.......750.........760.........770.........780.........790.........800.........810
.........820.........830.........840.........850.........860.........870........
.880.........890.........900.........910.........920.........930.........940......
...950.........960.........970.........980.........990........
999.

Done.
plot(F_CK.csr)

TAMPINES PLANNING AREA

Computing F-Function Estimation

F_tm = Fest(childcare_tm_ppp, correction = "best")
plot(F_tm)

Performing Complete Spatial Randomness Test

Ho = The distribution of childcare services in Tampines is randomly distributed.

H1= The distribution of childcare services in Tampines is not randomly distributed.

The null hypothesis will be rejected if the p-value is smaller than alpha value of 0.001.

F_tm.csr <- envelope(childcare_tm_ppp, Fest, correction = "all", nsim = 999)
Generating 999 simulations of CSR  ...
1, 2, 3, ......10.........20.........30.........40.........50.........60..
.......70.........80.........90.........100.........110.........120.........130
.........140.........150.........160.........170.........180.........190........
.200.........210.........220.........230.........240.........250.........260......
...270.........280.........290.........300.........310.........320.........330....
.....340.........350.........360.........370.........380.........390.........400..
.......410.........420.........430.........440.........450.........460.........470
.........480.........490.........500.........510.........520.........530........
.540.........550.........560.........570.........580.........590.........600......
...610.........620.........630.........640.........650.........660.........670....
.....680.........690.........700.........710.........720.........730.........740..
.......750.........760.........770.........780.........790.........800.........810
.........820.........830.........840.........850.........860.........870........
.880.........890.........900.........910.........920.........930.........940......
...950.........960.........970.........980.........990........
999.

Done.
plot(F_tm.csr)

5.3 Analyzing Spatial Point Process Using K-Function

The K-Function measures the number of events found up to a given distance of any particular event.

In this section, we are using:

  • Kest() of the spatstat package to compute the G-Function Estimation

  • envelope() of the spatstat package to perform Monte Carlo Simulation Tests.

CHOA CHU KANG PLANNING AREA

Computing K-Function Estimation

Observe that we are following the same process as the previous sections, just with a different estimation function.

K_ck = Kest(childcare_ck_ppp, correction = "Ripley")
plot(K_ck, . -r ~ r, ylab= "K(d)-r", xlab = "d(m)")

Performing Complete Spatial Randomness Test

We will conduct a hypothesis test to confirm the observed spatial patterns above.

Ho = The distribution of childcare services in Choa Chu Kang is randomly distributed.

H1= The distribution of childcare services in Choa Chu Kang is not randomly distributed.

The null hypothesis will be rejected if the p-value is smaller than alpha value of 0.001.

The code chunk below is used to perform the hypothesis testing.

K_ck.csr <- envelope(childcare_ck_ppp, Kest, nsim = 99, rank = 1, glocal=TRUE)
Generating 99 simulations of CSR  ...
1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20,
21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60,
61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80,
81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 
99.

Done.
plot(K_ck.csr, . - r ~ r, xlab="d", ylab="K(d)-r")

TAMPINES PLANNING AREA

Computing K-Function Estimation

K_tm = Kest(childcare_tm_ppp, correction = "Ripley")
plot(K_tm, . -r ~ r, 
     ylab= "K(d)-r", xlab = "d(m)", 
     xlim=c(0,1000))

Performing Complete Spatial Randomness Test

We will conduct a hypothesis test to confirm the observed spatial patterns above.

Ho = The distribution of childcare services in Tampines is randomly distributed.

H1= The distribution of childcare services in Tampines is not randomly distributed.

The null hypothesis will be rejected if the p-value is smaller than alpha value of 0.001.

The code chunk below is used to perform the hypothesis testing.

K_tm.csr <- envelope(childcare_tm_ppp, Kest, nsim = 99, rank = 1, glocal=TRUE)
Generating 99 simulations of CSR  ...
1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20,
21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60,
61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80,
81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 
99.

Done.
plot(K_tm.csr, . - r ~ r, 
     xlab="d", ylab="K(d)-r", xlim=c(0,500))

5.4 Analyzing Spatial Point Process Using L-Function

When I try to research on what exactly the purpose of L-Function is, it will usually just say that even mathematicians cannot come up with an exact definition of the function.

In this section, we are using:

  • Lest() of the spatstat package to compute the G-Function Estimation

  • envelope() of the spatstat package to perform Monte Carlo Simulation Tests.

CHOA CHU KANG PLANNING AREA

Computing L-Function Estimation

L_ck = Lest(childcare_ck_ppp, correction = "Ripley")
plot(L_ck, . -r ~ r, 
     ylab= "L(d)-r", xlab = "d(m)")

Performing Complete Spatial Randomness Test

We will conduct a hypothesis test to confirm the observed spatial patterns above.

Ho = The distribution of childcare services in Choa Chu Kang is randomly distributed.

H1= The distribution of childcare services in Choa Chu Kang is not randomly distributed.

The null hypothesis will be rejected if the p-value is smaller than alpha value of 0.001.

Monte Carlo Test:

L_ck.csr <- envelope(childcare_ck_ppp, Lest, nsim = 99, rank = 1, glocal=TRUE)
Generating 99 simulations of CSR  ...
1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20,
21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60,
61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80,
81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 
99.

Done.
plot(L_ck.csr, . - r ~ r, xlab="d", ylab="L(d)-r")

TAMPINES PLANNING AREA

Computing L-Function Estimation

L_tm = Lest(childcare_tm_ppp, correction = "Ripley")
plot(L_tm, . -r ~ r, 
     ylab= "L(d)-r", xlab = "d(m)", 
     xlim=c(0,1000))

Performing Complete Spatial Randomness Test

We will conduct a hypothesis test to confirm the observed spatial patterns above.

Ho = The distribution of childcare services in Tampines is randomly distributed.

H1= The distribution of childcare services in Tampines is not randomly distributed.

The null hypothesis will be rejected if the p-value is smaller than alpha value of 0.001.

Monte Carlo Test:

L_tm.csr <- envelope(childcare_tm_ppp, Lest, nsim = 99, rank = 1, glocal=TRUE)
Generating 99 simulations of CSR  ...
1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20,
21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60,
61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80,
81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 
99.

Done.

Then, plot the model output by using the code chunk below.

plot(L_tm.csr, . - r ~ r, 
     xlab="d", ylab="L(d)-r", xlim=c(0,500))